Sloppiness, robustness and evolvability in systems biology 
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The functioning of many biochemical networks is often robust - remarkably stable under changes 
in external conditions and internal reaction parameters. Much recent work on robustness and evolv- 
ability has focused on the structure of neutral spaces, in which system behavior remains invariant to 
mutations. Recently we have shown that the collective behavior of multiparameter models is most 
often sloppy: insensitive to changes except along a few 'stiff' combinations of parameters, with an 
enormous sloppy neutral subspace. Robustness is often assumed to be an emergent evolved property, 
but the sloppiness natural to biochemical networks offers an alternative non-adaptive explanation. 
Conversely, ideas developed to study evolvability in robust systems can be usefully extended to 
characterize sloppy systems. 



Introduction 

Robustness and evolvability are major themes of sys- 
tems biology, have been the subject of several recent 
books and reviews [H El El SI [5] , and have been discussed 
alongside related phenomena such as canalization, home- 
ostasis, stability, redundancy, and plasticity [6l [71 El E] . 
Broadly construed, "robustness is the persistence of an 
organismal trait under perturbations" , which requires 
the specification of both traits of interest and pertur- 
bations under consideration. Recent work in systems 
biology has sought to distinguish between environmen- 
tal robustness (e.g., temperature compensation in cir- 
cadian rhythms [lQiilHll2j) and mutational robustness 
(e.g., parameter insensitivity in segment polarity pattern- 
ing [I3l[l4]). Mutational robustness has a subtle relation 
to evolvability; while allowing survival under genetic al- 
terations, robustness might seem to reduce the capac- 
ity for evolutionary adaptation on multigeneration time 
scales HE]. 

Earlier robustness work focused on feedback and con- 
trol mechanisms [151 [161 |l71 [I8l [191 [20]. Much re- 
cent work emphasizes neutral spaces and neutral net- 
works: large regions in the space of sequences, param- 
eters, or system topologies that give rise to equivalent 
(or nearly equivalent) phenotypic behaviors. Neutral 
spaces have been explored most extensively in the con- 
text of RNA secondary structure, where large neutral 
networks of RNA sequences (genotypes) fold into iden- 
tical secondary structures (phenotypes) [SI [2ll [22l [23] . 
More recently, similar ideas have been applied to neu- 
tral spaces underlying the robustness of gene regulatory 
networks [24l [25l [26] , where different network topologies 
(genotypes) can result in identical gene expression pat- 
terns (phenotypes). Nontrivial niches in sequence spaces 
are also seen to emerge in molecular discrimination, a 
problem where neutral networks allow for biological com- 



munication in the presence of uncertainty akin to that 
found in engineered error-correcting codes [27 . Func- 
tional redundancies and degeneracies arise at many lev- 
els of biological organization [28], and it is an impor- 
tant open question as to how neutrality, redundancy, and 
robustness at different levels are organized and coupled 
across scales. 

Despite these advances in understanding neutral net- 
works connecting genotypes in discrete spaces (e.g., se- 
quences), much of systems biology is focused on chemical 
kinetic networks that are parameterized by continuous 
parameter spaces. Often one is interested in the steady- 
state behavior of a dynamical system, or in the input- 
output response relating only a subset of the chemical 
species of a network. In principle, however, one must 
characterize the full dynamical behavior of a network, in 
part because any given network may be coupled in un- 
known ways to other subsystems that are not included 
in the model. To more clearly delineate distinct levels of 
biological organization, we have chosen to refer the space 
of continuous kinetic parameters as a "chemotype" [29] . 
and to the full dynamical response of a system as its "dy- 
natype" (Figure [T]). The chemotype-to-dynatype maps 
of interest here are embedded within larger genotype-to- 
phenotype maps, with chemotypes emerging from lower- 
level processes, and dynatypes contributing to pheno- 
types and ultimately fitnesses on which selection acts. 
Recently, there has been increased interest in character- 
izing the parametric sensitivity of the dynamics of bio- 
chemical network models, for two important reasons: (1) 
to probe system robustness by quantifying the size and 
shape of chemotype spaces that leave system behavior 
unchanged, and (2) to characterize system behavior and 
uncertainties for which precise values for rate constants 
and other kinetic parameters are typically not known. 

Parameter estimation in multiparameter models has 
long been known to be ill-conditioned: the collective be- 
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havior usually cannot be used to infer the underlying con- 
stants. Recent work has shown that these models share 
striking universal features [30l ISTJ [32l |33] , a phenomenon 
that we have labeled "sloppiness" (see Figures [l]and[2|. 
Sloppiness refers to the highly anisotropic structure of pa- 
rameter space, wherein the behavior of models is highly 
sensitive to variation along a few 'stiff' directions (combi- 
nations of model parameters) and more or less insensitive 
to variation along a large number of 'sloppy' directions. A 
nonlinear least-squares cost function can be constructed: 

i * i 

where = {x{0) — Xi)/(Ji is the residual describing 
the deviation of a dynamical variable x from its mea- 
sured values Xi with uncertainty cr^. This cost reflects 
how well a model with a given set of parameters 
fits observed experimental data. Parametric sensitivi- 
ties of the model are encoded in the Jacobian matrix 
J = dri/dOj. The curvature of the cost surface about 
a best fit set of parameters is described by the Hessian 
Hmn = d'^C / dOjnOn (oi" l^s approximation, the Fisher 
Information Matrix J^J). Stiff and sloppy directions 
are conveniently measured using an analysis of eigen- 
values An of the Hessian H (Figure [3| ; large eigenval- 
ues correspond to stiff directions. For a broad range 
of multiparameter models (e.g., sixteen models drawn 
from the systems biology literature [32 and models from 
quantum Monte Carlo, radioactive decay, and polyno- 
mial fitting [34 ) these eigenvalues are roughly uniformly 
spread over many decades, with many sloppy directions 
a thousand times less well determined than the stiffest, 
best constrained parameter combinations. Two conse- 
quences are that useful model predictions can be made 
even in the face of huge remaining parameter uncertainty, 
and conversely that direct measurements of the param- 
eters can be inefficient in making more precise predic- 
tions [32 . Random matrix theory can be used to de- 
velop insight into the source of this type of eigenvalue 
spectrum and the nature of redundancies that appear to 
underly sloppiness [34 . Our open-source code SloppyCell 
(http://sloppycell.sourceforge.net) provides tools 
for exploring parameter space of systems biology mod- 
els [35]. 

Others have recently addressed similar questions mo- 
tivated by the lack of detailed information about kinetic 
parameters. These include: the inference of probabilis- 
tic statements about network dynamics from probability 
distributions on parameter values [36 ; the use of "struc- 
tural kinetic modeling" to parameterize the Jacobian ma- 
trix J and thereby probe ensembles of dynamical behav- 
iors [371 EH] ; the construction of convex parameter spaces 
("k-cones") containing all allowable combinations of ki- 
netic parameters for steady-state flux balance [39]; the 
use of ideas from control theory, worst-case analysis and 
hybrid optimization to measure the robustness of net- 
works to simultaneous parameter variation [40], and ex- 
ploration of correlated parameter uncertainties obtained 




FIG. 1: Sloppiness in the mapping of chemotypes to 
dynatypes. It is natural, at least for cellular regulation 
and metabolic networks, to refine the traditional dichotomy 
of genotype G to phenotype P by adding two intermediate 
levels of description, G ^ C ^ D ^ P. Here C is the 
chemotype 29^, a continuous description of the behavior in 
terms of chemical reaction parameters (reaction rates, bar- 
riers and prefactors, or Michaelis-Menten parameters). D is 
the dynatype, meant to describe the dynamical responses of 
the cell (usually the time series of all species in response to se- 
lected stimuli, often taken from experimental measurements). 
Mutations about a particular chemotype 6 occupy a region in 
chemotype space (here a circle of radius S), whose image in 
dynatype space is given by the local Jacobian J of the map- 
ping: mutations along stiff directions in chemotype space will 
yield large changes in dynatype, while mutations along sloppy 
directions will lead to small dynamical changes. Conversely, 
a population of individuals sharing nearly the same dynatype 
r (here a sphere of radius e) will occupy a distorted region in 
chemotype space, with large variations in reaction parameters 
possible along sloppy directions (gray ellipse). 



via global inversion [41]. 

Can we connect sloppiness to robustness and evolv- 
ability? It is our contention that sloppiness - the highly 
anisotropic structure of neutral variation in the space of 
chemotypes - has important implications for how one 
characterizes robustness in systems biology models. In 
addition, insights developed in the study of robustness 
and evolvability suggest new and potentially useful ways 
of analyzing and interpreting sloppiness. 



Environmental robustness and sloppiness 

Organisms must thrive under many environmental con- 
ditions: changing temperatures, salt concentrations, pH, 
nutrient densities, etc. Many organisms have explicit 
control mechanisms to keep their internal state insensi- 
tive to these external changes - these control mechanisms 
(homeostasis, adaptation, etc.) have been a historical fo- 
cus in the robustness literature [15] |43] . For variations 
in temperature, however, many organisms do not have 
such homeostatic control (with the exception of birds, 
mammals, and some plants) and must instead cope with 
the exponential Arrhenius temperature dependence of all 
their reaction rates by some sort of compensatory mech- 
anism [44 . 

The prototypical example of temperature compensa- 
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FIG. 2: Sloppy parameter distributions: dependence 
on external conditions. Shown is a two-dimensional view 
of the parameter sets (free energy barriers and prefactors) 
that accurately predict the experimental phosphorylation dy- 
namics 11 in a 36-parameter subnetwork of a model of cir- 
cadian rhythms 12 , within a harmonic approximation (see 
Supplemental Material) . Shown are parameters valid at three 
different temperatures (colors) and valid for all temperatures 
simultaneously (black). The plot shows one 'stiff' direction 
in parameter space for each temperature which is tightly con- 
strained by the data, and one 'sloppy' direction which has 
relatively large variations without change in behavior. Most 
of the 34 other directions in parameter space not shown are 
sloppy; the two-dimensional view was chosen to best align 
with the stiffest direction for each of the four ensembles. The 
black region models organisms that are robust to temperature 
changes in this range. The acceptable region rotates and shifts 
with temperature, but the sloppiness allows different temper- 
atures to intersect (robust temperature compensation) even 
though all rates are strongly temperature dependent. 



tion is the 24-hour period of circadian rhythms [10]. Re- 
cent experiments have succeeded in replicating the cir- 
cadian control network of cyanobacteria in the test tube 
using three Kai proteins, whose degree of phosphoryla- 
tion oscillates with a temperature-compensated period in 
the range of 25 to 35° C. In addition, the phosphorylation 
dynamics of KaiC alone is found to be unchanged as the 
temperature varies in the same range [llj . This has been 
cited as a plausible explanation for the observed temper- 
ature compensation in the full network, presuming that 
all other rates are fast [12 and hence irrelevant to the 
period. (At least one other explanation of temperature 
compensation [45 also relies on constraining most rates 
to be irrelevant). Narrowing our focus to the KaiC phos- 
phorylation subnetwork, however, still leaves the non- 
trivial task of explaining its temperature compensation 
mechanism, since estimated energy barriers [46 suggest 
that phosphorylation rates should be twice as fast at the 
higher temperature. 

The dynamics of KaiC phosphorylation have been 
modeled using six phosphorylation sites and two confor- 
mational states (active and inactive) [12 . If each of the 
18 rates in this model roughly double between 25 and 
35° C, can we adjust the corresponding energy barriers 
and prefactors such that the resulting net phosphoryla- 
tion dynamics is temperature-independent? 

Figure [2] shows a two-dimensional view of the accept- 
able parameter sets in the resulting 36-dimensional space 
of energy barriers and prefactors, explored in the har- 
monic approximation (see Supplemental Material). No- 
tice that the region of acceptable parameters rotates and 




FIG. 3: Sloppy model eigenvalues. Shown are the eigen- 
values of the approximate Hessian J^J for the goodness-of- 
fit C(d) (Equation [T]) about the best fit. Large eigenvalues 
correspond to stiff directions; others are sloppy. Notice the 
enormous range on this logarithmic scale; not all eigenvalues 
(ranging down to 10~^°) are depicted. 

• Columns KaiC 30 and KaiC All are for the KaiC phosphory- 
lation dynamics model (Figure[3|, showing T — 30° C (yellow 
region in Figure [2]) and simultaneous fits for all temperatures 
(black region). Notice that the 'robust' simultaneous fit has 
roughly one more stiff direction than the single temperatures. 

• The SP and SP PC A columns are for the segment polarity 
model [I3l[42]. SP is an eigenvalue analysis about one of the 
acceptable parameter sets, showing parameters that keep the 
behavior (dynatype) of the entire network preserved (time se- 
ries for all components under all experimental conditions) . SP 
PC A is a principal components analysis of the segment po- 
larity ensemble that yields the wild- type phenotype, with pa- 
rameters restricted to a relatively small range (roughly three 
decades each). Most directions in SP are sloppy enough to 
have fluctuations larger than the sampled phenotype box in 
SP PC A] the sloppy dynatype SP already explains the ro- 
bustness to all but a few stiff directions in parameter space. 
Conversely, the sensitivity of the dynatype SP to a few stiff 
directions does not preclude phenotypic robustness in those 
directions for SP PC A] the dynatype (all dynamical evolu- 
tion) is far more restrictive than the phenotype (output pat- 
terning) . 

• PC 12 is for the EGF/NGF growth-factor signahng net- 
work [31] [32]; note that it too is sloppy. See Figure |4]for 
an analysis of evolvability and robustness for this model. 



shifts as the temperature changes. Notice also that the 
system is sloppy: Figure |2] shows one stiff direction that is 
highly constrained by the data and one sloppy direction 
that is largely unconstrained. The eigenvalue analysis 
in Figure |3] confirms that most directions in parameter 
space are sloppy and unconstrained. This provides a nat- 
ural explanation for robustness: the intersection of these 
large, flat hypersurfaces yields parameters that work at 
all temperatures. [32 In general, each external condition 
provides one constraint per stiff direction; since there are 
only a few stiff directions and many parameters in sloppy 
models, robust behavior under varying external condi- 
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tions is easily arranged. Indeed, Figure [3] shows that the 
robust, temperature-independent fits for the KaiC model 
are themselves a sloppy system. 



Chemotype robustness and sloppiness 

In addition to robustness to environmental perturba- 
tion, biological networks are often robust to mutational 
perturbations; they maintain their function in the face 
of mutations that change one or perhaps more of their 
underlying rate parameters, and thus change their loca- 
tion in chemotype space. Some authors have used this 
as a criterion for judging model plausibility [47|. The 
quintessential example of a system that is chemotypi- 
cally robust is the Drosophila segment polarity gene net- 
work. Early in development, this network generates a 
periodic macroscopic phenotype: a pattern of gene ex- 
pression across several cells that persists throughout de- 
velopment and guides later stages. Multiparameter mod- 
els of this network [131 ttH 113 SH] find that a surprisingly 
large fraction of randomly chose parameter sets gener- 
ate a pattern consistent with the observed patterning of 
three genes - the system exhibits chemotype robustness. 

In the context of sloppy models, we may define chemo- 
type robustness as the fraction of a given volume in pa- 
rameter/chemotype space C that maps into a functional 
region of behavior /dynatype space D (Figure [l]). This 
latter functional region represents behavior close to op- 
timum (or close to that measured experimentally). For 
simplicity, let us consider it to be a hypersphere of ra- 
dius e (i.e., a cost C{0) = < in Equation [l]) ; 
larger changes in behavior are considered significantly 
different, perhaps lowering the organism's fitness. The 
given volume in chemotype space C might be (as for 
the segment polarity network) a hypercube of parame- 
ter ranges deemed reasonable, or (as a simple model of 
mutations) a hypersphere; let its scale be given by 5. Our 
robustness is therefore the fraction of all points in the S- 
ball in C that map into the e-ball in D - in Figure [l] the 
fraction of the circle whose interior is colored gray. This 
fraction can be calculated (see Supplemental Material) 
and is approximately given by 

where Xcrit = This formula can be motivated by 

considering the robust subregion (gray needle intersect- 
ing the circle) to be a slab, with thickness e^/A^ along the 
eigendirection corresponding to each eigenvalue Ar,.[5Q] 
For sloppy directions with < = Xcrit ^ the slab 

is thicker than the circle and does not reduce the ro- 
bust fraction; for each stiff direction with A^ > Xcrit ^ the 
fractional volume is reduced roughly by a factor of the 
slab thickness over the sphere width ^, leading to 

Equation 



In their model of segment polarity, von Dassow et al. 
found that approximately one in 200 randomly chosen pa- 
rameter sets generated a wild-type expression pattern for 
three key genes [TT . This would naively seem amazing for 
a 48 parameter model like theirs; in an isotropic approxi- 
mation, each parameter would be allowed only 6% chance 
of changing the wild- type pattern (since 0.94^^ ~ 1/200). 
However, we have previously shown that the segment po- 
larity model is sloppy [32 . That is, going far beyond 
restricting the output phenotype, the dynamical evolu- 
tion of every component of the network is approximately 
preserved even with huge changes in parameter values: 
only a few stiff directions in chemotype space are needed 
to maintain the dynatype (see column SP in Figure [3| . 
Sloppiness hence provides a natural explanation for the 
wide variations in all but a few directions in parameter 
space. 

The success rate of one in 200 is not nearly as strik- 
ing if the dynamics is already known to be insensitive to 
all but perhaps four or five combinations of parameters: 
0.35^ X 1^3 - 1/200. Column SP PC A in Figure |3]fieshes 
this picture out with a principal components analysis 
(PGA) of the robust region seen in von Dassow et al.'s 
original model, reconstructed using Ingeneue [42 . Note 
that these PGA eigenvalues are cut off from below by the 
parameter ranges chosen by the original authors for ex- 
ploration (typically three decades per parameter) . While 
the overall scale of the dynatype sloppy-model eigenval- 
ues in SP and the phenotype eigenvalues in SP PC A can- 
not be directly compared, it is clear that the vast major- 
ity of sloppy-model eigenvalues are too small to constrain 
the parameters within the explored region. The model is 
robust in these directions not because of evolution and 
fitness, but because the dynamics of chemical reaction 
networks is mathematically naturally dependent only on 
a few combinations of reaction parameters. 



Robustness, evolvability, and sloppiness 

Mutational robustness of systems would seem to be 
at odds with an ability to adapt and evolve, since ro- 
bustness implies persistence of phenotype or function, 
which may inhibit the capacity for evolutionary change. 
The concept of neutral spaces has been used - most no- 
tably by Wagner and collaborators - to suggest a res- 
olution of this apparent paradox, as demonstrated in 
model systems exploring various genotype-to-phenotype 
maps [8[ [23[ [Ml ES] • The important insight is that neu- 
tral spaces and neutral networks enable systems to drift 
robustly in genotype space (i.e., without significant phe- 
notypic change), while encountering new and different 
phenotypes at various points along that neutral space. 
This insight results from a distinction between the ro- 
bustness and evolvability of any given genotype, and the 
robustness and evolvability of all genotypes consistent 
with a given phenotype [8 . 

Evolvability is postulated to refiect the range of possi- 



5 



by: 



Population 




ec(F,6f) = VF^JJ^FS 



(3) 



10 



10 10 10 

Evolvability e(F) 



10 



FIG. 4: Evolvability and robustness in a sloppy sys- 
tem. Evolvability distributions, and evolvability versus ro- 
bustness, for an ensemble of parameters for a model of an 
EGF/NGF signaling pathway fitted to experimental data in 
PC12 cells [32]. The histogram on the left is the distribu- 
tion of individual/chemotype evolvabilities ec{F,6a) (Equa- 
tion [3]) , as F (an evolutionary pressure in dynatype space) 
is randomly chosen in direction with uniform magnitude and 
Ocx varies over the ensemble. The histogram on the right is 
the corresponding distribution of population/dynatype evolv- 
abilities ed{F) (Equation [4| . Note that the population evolv- 
abilities are significantly higher than the individual ones. The 
inset plots the RMS individual chemotype evolvability EdOa) 
versus the robustness Rc{9oi) (Equation |2]) for the ensemble. 
{Xcrit is chosen as the fourth- stiff est eigenvalue at the best 
fit: see Supplemental Material). Note that, for each individ- 
ual, more robustness leads to less evolvability - individuals 
which rarely mutate to new forms can't evolve as readily. This 
need not apply to the population, insofar as we expect robust 
dynatypes to explore larger regions of parameter/chemotype 
space, and thus the ratio of dynatype to chemotype evolvabil- 
ity to increase with increasing robustness. 



ble different phenotypes that are possible under geno- 
typic mutation. How does the sloppy connection be- 
tween parameters and behavior impinge on the question 
of evolvability? Translating previous work on discrete 
genotype and phenotype spaces to the continuous spaces 
of chemotypes and dynatypes is nontrivial. Since the 
dimensionality of the space of chemotypes is less than 
that of dynatypes, the volume of dynatype space acces- 
sible under changes in chemotype is zero, i.e., lies on a 
lower-dimensional subspace. To develop a sensible defi- 
nition of evolvability in such systems, we postulate forces 
F in dynatype space (Figure [T]) that reflect evolutionary 
pressures due to changes in the environment, such that 
a change r in dynatype leads to a change r • F in fitness. 
An organism's evolvability is related to its capacity to re- 
spond to external forces through appropriate mutations 
in chemotype. 

For a given force F, the maximum fitness change 
among mutations of size d in chemotype space is given 



which we call the chemotype evolvability distribution (see 
Supplemental Material). Refs. [31] and [32] generate en- 
sembles of parameters (chemotypes) consistent with a 
given dynatype for an EGF/NGF signaling pathway in 
PC 12 cells, where the dynatype is constrained to fit avail- 
able experimental data. (The PC 12 network is sloppy, see 
Figure 3]) Each member of such an ensemble Ooc has a 
Jacobian Jq,. As in Ref. [8 , which distinguishes between 
genotype and phenotype evolvability, we can distinguish 
between the chemotype eciF^Oa) and dynatype 



erf(F) = maxec(F,^c 



(4) 



evolvability distributions. The first gives the distribution 
of adaptive responses to F of individual chemotypes in a 
population, while the second gives the optimal response 
within the population. Figure [4] shows the chemotype 
and dynatype evolvability distributions, generated using 
the PC 12 ensemble of Ref. [32 and a uniform distribu- 
tion of force directions F in dynatype space. Within 
a population sharing the same behavior, we find sub- 
stantial variation of accessible behavior changes, leading 
to a substantially larger population (dynatype) evolv- 
ability than individual (chemotype) evolvability. This 
echoes the finding of Wagner that phenotype evolvability 
is greater than genotype evolvability for RNA secondary 
structures [8^. 

It is natural to define an overall evolvability as the 
root-mean-square average of the evolvability distribution 
over a spherical distribution of environmental forces F in 
dynatype space: 



(5) 



and correspondingly for the overall RMS dynatype evolv- 
ability. The inset to Figure |4] shows that the chemo- 
type evolvability decreases as the chemotype robustness 
increases, closely analogous to Wagner's discovery that 
genotype evolvability decreases as genotype robustness 
increases, except that his plot averages over phenotypes 
while ours represents variation within a dynatype. Thus 
we reproduce Wagner's observation [8] that individual 
evolvability decreases with robustness and that popu- 
lation evolvability is significantly larger than individual 
evolvability. [51 



Conclusion 

Our previous work aimed at developing predictive sys- 
tems biology models in the face of parametric uncertainty 
has led us to formulate a theory of sloppiness in multi- 
parameter models. The picture that emerges from this 
theory is of a highly anisotropic neutral space in which 
variation in parameters (chemotypes) can leave system 
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behavior (dynatypes) unchanged. This picture is rem- 
iniscent in many ways to the notion of neutral spaces 
and neutral networks that has been developed to explore 
the robustness and evolvability of biological systems. We 
have been motivated by those ideas to here reconsider 
sloppiness within that context, both to highlight implica- 
tions of sloppiness for the study of robustness and evolv- 
ability, and to identify new methods for analyzing sloppy 
systems. 
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Contents 

The contents of the supplemental material are orga- 
nized corresponding to the order of the main text. In- 
cluded in the supplemental material are derivations of 
mathematical results and details of the specific models 
mentioned in the main text. 

We have also posted the data files and computer 
codes for the models discussed, at http: //www. lassp. 
Icornell . edu/ sethna/Sloppy , For the KaiC, PC12, and 
segment polarity models, this includes: 

1. Equations in I^TfeX, Python, and C 

2. SBML (system biology markup language) files 

3. Parameter ensembles 

4. Best-fit Hessian and J-^J, and their eigenvectors 
and eigenvalues 



Introduction 

Hessian at best fit parameters 

In the introduction we mention that "the curvature 
of the cost surface about a best fit set of parameters is 
described by the Hessian Hmn-^^ Examining the behavior 
of Hmn is a standard method for nonlinear least squares 
models when fitting data. Formally, Hmn is written as: 



Hrn/n. — 



= E 



drj drj 



(SI) 



If the model fits the data well so that r 
Ref. [32 in the main text) then 

dri dr, 



(or perfectly. 



{J' J). 



(S2) 



If H and the cost are used (as in this review) to describe 
changes in model behavior from then r = at ^* and 
Equation (S2) is exact. Notice also that H refiects the 



sensitivity of the fit to changes in parameters; in fact, 
its inverse is the covariance matrix. The diagonal ele- 
ments of the covariance matrix are proportional to the 
uncertainties in the parameters, while the off-diagonal 
elements are estimates of parameter uncertainty correla- 
tions. 



Figure [7|' Sloppiness in the mapping of chemotypes to 
dynatypes 

Shown in Figure [l] of the main text is the mapping 
of chemotypes C to dynatypes D. The mapping be- 
tween C and D is described with J and "J~^". It is 
typical that dim(C) <C dim(D), since there are typically 
more data points constraining the dynatype than there 
are parameters defining a chemotype. Therefore, the in- 
verse of J is not well-defined. In Figure [l] the gray el- 
lipse in C represents the inverse image of the e-ball, B^, 
in D under J. That is, "J~^" acting on is the set 
{ceC s.t. J -ce Be}. 

Note also that the stiff and sloppy eigendirections in 
C and their images in D can be described by the sin- 
gular value decomposition of the Jacobian J. Since 
are eigenvalues of J^J, V%i are the singular values of 
J. Furthermore, writing J = U^V^^ we see that 
the columns of V are stiff/sloppy eigenparameters in C 
(shown in red in the figure), and the columns of U are 
images of stiff and sloppy eigenparameters (divided by 
An) in D. 



Environmental robustness and sloppiness 

Figure ^ Sloppy parameter distributions: dependence 
on external conditions 

In Figure |2] of the main text, the plane onto which the 
ensembles are projected is the one that aligns best with 
the stiffest eigenparameter of each of the four ensembles. 
To accomplish this, the vertical and horizontal axes in 
Figure |2] are, respectively, the first and second singular 
vectors in the singular value decomposition of the set 
of stiffest eigenparameters {vq^, Vq^, Vq^, Vq^^^}. In a way 
analogous to principal components analysis, this gives the 
plane that passes through the origin and comes closest to 
passing through the heads of unit vectors pointing in the 
stiffest eigendirections. 

Each ensemble of parameter sets shown in Figure [2] is 
chosen from the probability distribution corresponding 
to the local quadratic approximation of the cost near the 
best-fit parameters ^*: 



P(6/* + AO) (X exp{-A0J^JA0/2), 



(S3) 



This local approximation to the cost was used to generate 
the ensembles instead of the full nonlinear cost function 
due to difficulties in generating equilibrated ensembles: 
the thin curving manifolds of allowable chemotypes for 
sloppy models can be notoriously difficult to populate. 
But this is not impossible; efforts are still underway, and 
if equilibrated ensembles are found, they will be posted 
to the website mentioned above. 

KaiC phosphorylation subnetwork model 

In the main text, we use as an example a portion of 
the circadian rhythm model presented in Ref. [12^ of the 
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FIG. SI: KaiC phosphorylation subnetwork. This 
schematic depicts the KaiC network used as an example in 
the main text. It is a portion of the full circadian rhythm 
model presented in Ref. [12] of the main text. The num- 
bers represent the degree of phosphorylation, and the two 
columns represent two different conformational states, "ac- 
tive" and "inactive." The labels on the arrows represent 
reaction rates for changing among the phosphorylation and 
conformation states. Each conformation state has one phos- 
phorylation and one dephosphorylation rate, independent of 
the degree of phosphorylation. Each of the 14 "flip" rates 
between conformational states (bi and fi) is allowed to vary 
independently. This gives a total of 18 reaction rates. 



FIG. S2: KaiC phosphorylation network: 

temperature-compensated output. Shown is the 
net phosphorylation of KaiC over time, comparing experi- 
mental data (circles with error bars, from Ref. [11 in the 
main text) with output from an ensemble of chemotypes 
(filled colored regions, showing the mean plus or minus one 
standard deviation over the ensemble for the net phospho- 
rylation at each time-point). Different colors correspond to 
different temperatures: blue = 25°, green = 30°, red = 35°. 
Note that the chemotypes describe the data well at all three 
temperatures, even though the rates are strongly dependent 
on temperature. 



this prior as a quart ic in log E: 



main text. We implement the subnetwork that van Zon 
et al. hypothesize must have intrinsically temperature- 
independent rates: that which controls the phosphoryla- 
tion of KaiC alone. This subnetwork models the experi- 
mental measurements of KaiC phosphorylation in the ab- 
sence of KaiA and KaiB (Ref. [11 of the main text), in 
which the phosphorylation of KaiC does not oscillate, but 
decays at a temperature-compensated rate in the range 
from 25 to 35° C (see circles in Figure S2). 

The subnetwork involves an active and inactive state 
of KaiC, along with six phosphorylation sites for each 
state, as depicted in Figure [Si] Including forward and 
backward "flip" rates between active and inactive states 
along with (de) phosphorylation rates that are each con- 
stant for the two states, there are 18 independent rates. 
To assess the temperature dependence, we assume that 
each transition rate follows an Arrhenius law, with con- 
stant energy barrier E and prefactor a: the ith rate is 
^,^Ei/kT ^ This then gives a 36-dimensional chemotype 
space in which to search for solutions. 

Temperature-independent solutions can be trivially 
found in this space if the energy barriers are chosen to 
be small, since this produces rates that are inherently 
weakly dependent on temperature. In order to avoid 
this trivial temperature compensation, we apply a prior 
that favors solutions with phosphorylation energy barri- 
ers near the expected Eq = 23 /cT, similar to those found 
in other kinases (Ref. [46 in the main text) and appropri- 
ate for reactions that break covalent bonds. We choose 



25 

y 



log 



(S4) 



The form was chosen to severely penalize barriers less 
than 10 /cT, but to be reasonably flat around Eq] other 
prior choices would presumably perform similarly. 

Using this method, we find that it is possible to fit the 
experimental data even with (de)phosphorylation rates 
that are strongly temperature-dependent. The phos- 
phorylation and dephosphorylation rates that provided 
a best fit to all temperatures simultaneously were all 
above 21 kT. We used Bayesian Monte- Carlo sampling of 
chemotype space to create an ensemble of parameter sets 
that each produce phosphorylation dynamics that match 
the experimental data at 25, 30, and 35° C. As explained 
above, our ensemble has not yet sampled all the space 
available, but we still find many such acceptable chemo- 
types. The minimum (de)phosphorylation rate for the 
ensemble was just under 10 /cT, so the prior worked as 
designed to confine the barriers to physically reasonable 



values. Figure S2 shows the output of the model over 
this ensemble of parameter sets compared with the ex- 
perimental data from Ref. [11^ of the main text. 

We mention in a footnote that, in our model, "suc- 
cessful chemotypes favor dephosphorylation in the active 
state and phosphorylation in the inactive state." This 
can be seen in the ratio of phosphorylation to dephos- 
phorylation rates, shown in Figure [S3j for the ensemble 
of successful chemotypes. Note that most members of the 
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FIG. S3: KaiC phosphorylation network: 

temperature-compensation mechanism. This plot 
shows the ratios of phosphorylation rates to dephosphoryla- 
tion rates for the active and inactive states - the distribution 
of kps/kdps is shown in blue for the active state, and the dis- 
tribution of k^s/kdps is shown in green for the inactive state 



(see Figure SI for definitions of rate constants). The distribu- 
tion is over the same (non-equilibrated) ensemble as was used 
to generate Figure |S2] Note that phosphorylation is favored 
in the inactive state, while dephosphorylation is favored in 
the active state. This suggests a temperature-compensation 
mechanism, as described in the text. 



ensemble have an inactive state with higher phosphoryla- 
tion rate than dephosphorylation, and vice versa for the 
active state. This matches with an intuitive temperature- 
compensation mechanism: with flip rates that are also 
temperature-dependent, higher temperatures can lead to 
more KaiC being in the inactive state, leading to a slower 
overall decay in phosphorylation that compensates for 
the speedup in reaction rates. 

Figure\9^ Sloppy model eigenvalues 

The PCA shown in Figure |3] column SP PCA was pro- 
duced after taking logarithms of the parameter values 



that von Dassow et al. used in their analysis. This 
measures parameter fluctuations in terms of fractional 
changes in parameter, rather than absolute sizes of fluc- 
tuations - allowing fluctuations in parameters with differ- 
ent units, for example, to be compared. The parameters 
used in column SP were chosen (logarithmic or other- 
wise) as defined by the original authors. Taking loga- 
rithms and/or changing units does not typically change 
the qualitative spectra of sloppy models, as their spectra 
already span so many decades. 



Chemotype robustness and sloppiness 

Derivation of robustness equation 

In the main text (MT), the robustness is defined as 



Rn 



n 



A^>Ac 



(MT[2| 



We now proceed to derive this result. We measure ro- 
bustness as the fraction of mutations of a given size 5 
in C (chemotype space) that do not change the behav- 
ior beyond a given threshold (survival after a mutation), 
which we designate as an e-ball around the optimum in 
D (dynatype space). Therefore we want an estimate of 
the fraction of the (5-ball in C that maps into the e-ball 
m D. It is difficult to calculate this geometrically, since 
we would need to find the volume of an ellipsoid inter- 
secting a sphere. Fortunately, for sloppy systems, the 
\i vary over many orders of magnitude, so we can sim- 
plify the calculation by smearing the (5-ball and e-ball 
into Gaussians. Namely, we say a mutation in C has 
probability e'^^^^^ 1'^^^ / , and the probability of 

2 In 2 

"survival" in D is given by ' . We then measure the 
robustness as the overall probability P(J, e) of surviving 
after a mutation: 



Rc = P{S,e) 

^\ f dAO exp(-(A6/)V2^^)exp(-(A6/)^J^J(A6/)/2e2 
27tS J Jc 



TT , ^ =. (S5) 



For sloppy systems, A varies over many orders of magni- Therefore, using our definition Xcrit = e^/J^ we can ap- 
tude. Notice that if <C its component in the 

product will be close to 1, and if A^ ^ ^ 1^^ ^ we can ap- 
proximate the components in the product as yje^ j^^X^ 
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proximate this formula as: 



RMS dynatype evolvability 
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^crit 
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with smah corrections for eigenvalues ^ e'^/S'^. Since 
this result agrees with the "slab" argument given in the 
main text for hard walls, we see that hard e-balls and 
hard J-balls will have approximately the same amount of 
overlap as Gaussians. 



Robustness, evolvability, and sloppiness 

Derivation of chemotype evolvability 

In the main text, we provide a formula for the "maxi- 
mum fitness change among mutations of size S in chemo- 
type space" 



ec(F,6/) = Vf^JJ^FS 



(MT[3| 



which we derive here using a Lagrange multiplier. To 
derive this, we use the definition of the chemotype evolv- 
ability as the maximum response r • F in for moves in 
C of size \ Ae\ = 6: 



eJF.O) = max (r • F). 



(S7) 



Next, notice that 

r-F= (JA0)-F = ^^FiJi„A^„. (S8) 

i a 

We find the optimal using a Lagrange multiplier A. 
With {AO)'^ = S'^ as our constraint, we maximize 

F,j,^Aea^A{{Aef-s^) = F,j,^Aea^A{AepAef3-s^) 

(S9) 

where we use the Einstein summation convention (sum- 
ming over repeated indices). Differentiating with respect 
to AOa^ we can find the change A9^^^ giving the maxi- 
mum response: 



and hence 
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\nmax _ ^ J^JQ 



which implies 
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4<52 



Therefore, the evolvabihty is: 

ee(F,0) = F,Ji„AC°" 
F^JJ^F 



p. J. p. T. 

~ 2A 
5 = ^/YTjjTY5. 
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In Equation ([5|, to measure overall evolvability, we 
defined Ec{Oa) as a root-mean-square (RMS) average 
over a uniform (hyper) spherical distribution of environ- 
m ental forces F in dynatype space. We use the RMS 
^/{e^JF\Oa)^ rather than the average (edF^Oa)) be- 
cause the RMS definition has an elegant result in terms 
of the eigenvalues A^ of J: 

2^ fd^F 

i 

_ Tr{rj){F') ^ Tr{H){F') ^,^ ^^^^^ 



N 



N 



Therefore, the overall evolvability is directly related to 
the trace of the Hessian: 



Ec{Oa) 



Tr{H){F^) 
TV 



S. 
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Our measures of robustness and evolvability depend 
upon our level of description, just as for Wagner's geno- 
type and phenotype evolvabilities of RNA sequences 
(Ref. [8 of the main text). Our choice of an isotropic dis- 
tribution of selective dynatype forces F is not intended 
as an accurate representation of actual selective forces at 
the phenotype level, but as an exhaustive study of all 
possible forces at the dynatype level of description. 

Information about phenotypic selective pressures 
might suggest a different distribution of dynatype forces 
F. Indeed, this formalism provides a mechanism for cou- 
pling maps across scales, which is an important unsolved 
problem. Just as the genotype-to-chemotype (G C) 
and chemotype-to-dynatype {C D) maps are many-to- 
one, so is the dynatype-to-phenotype map {D ^ P). In 
the segment polarity model, for example, one might con- 
strue the phenotype as the steady-state pattern, whereas 
the dynatype will include information about all transient 
paths to that steady state. This is also closely analogous 
to measuring evolvability of RNA sequences by counting 
distinct folded structures (Ref. [8 of the main text), as 
many different structures may be equally nonfunctional 
at the higher level of biological phenotype. Ultimately, 
understanding the nature of the complex D ^ P maps 
will be required to estimate evolvability using more real- 
istic distributions of selective dynatypic forces F. 

Figure^ Evolvability and robustness in a sloppy system 

When calculating the chemotype robustness Rc^ we 
have a choice to make for the value of Xcru (see Equa- 
tion [2| . This choice corresponds to setting the ratio of 
the size of acceptable changes in dynatype e to the typi- 



cal size of mutations 5 in chemotype space: Ac 
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Equivalently, Xcrit sets a cutoff between stiff and sloppy 
eigenvalues, since we assume that, in D space, the im- 
age of the (5-ball fully overlaps with the e-ball in sloppy 
directions (with eigenvalues below Xcrit) ^ and it extends 
far beyond the edge of the e-ball in stiff directions (with 
eigenvalues above Xcrit)- 

In calculating Rc for the inset of Figure [4] in the main 
text, we chose Xcrit as the fourth stiffest eigenvalue of 
J^J at the best fit parameters. This matches with 
the idea that there are only a few stiff directions that 
appreciably constrain parameters in chemotype space: 
the eigenvalues are spaced by roughly factors of three 
(Figure |3|, meaning mutations in sloppier directions in 
chemotype space quickly become irrelevant in dynatype 
space. The choice of Xcrit within a reasonable range (be- 
tween, say, the second stiffest and eighth stiffest eigen- 
value of J^J) does not qualitatively change the plot of 
evolvability vs. robustness. 



